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AN OPTIMIZATION-BASED ATOMISTIC-TO-CONTINUUM 
COUPLING METHOD* 

DEREK OLSONt, PAVEL B. BOCHEVt, MITCHELL LUSKINt, AND ALEXANDER V. 

SHAPEEVt 

Abstract. We present a new optimization-based method for atomistic-to-continuum (AtC) cou- 
pling. The main idea is to cast the coupUng of the atomistic and continuum models as a constrained 
optimization problem with virtual Dirichlet controls on the interfaces between the atomistic and con- 
\^i tinuum subdomains. The optimization objective is to minimize the error between the atomistic and 

^~^ continuum solutions on the overlap between the two subdomains, while the atomistic and continuum 

(_^ force balance equations provide the constraints. Splitting of the atomistic and continuum problems 

fNj instead of blending them and their subsequent use as constraints in the optimization problem dis- 

> , tinguishes our approach from the existing AtC formulations. We present and analyze the method in 

O , the context of a one-dimensional chain of atoms modeled using a linearized two-body next-nearest 

^^ neighbor interactions. 

l^ 1. Introduction. Atomistic-to-continuum (AtC) coupling methods aim to com- 

bine the efficiency of continuum models such as PDEs with the accuracy of the atom- 
istic models necessary to resolve local features such as cracks or dislocations that can 
affect the global material behavior. Specifically, suppose that an atomistic model gives 
an accurate description of the true material behavior in a domain 17, but that this 
model is prohibitively expensive to solve on the whole domain. The core of AtC for- 
mulations is to keep this model only where the fully atomistic description is required 
Cw to accurately represent local features, while utilizing a more efficient continuum model 

H in the rest of H,. Existing AtC methods differ chiefly by the manner in which these 

'~^ models are joined together, which is also the main challenge in the AtC formulation. 

I To explain the main ideas we can consider a scenario where il is subdivided into 

^ an atomistic and a continuum subdomain, ila and flc, such that fla U Qc ~ ^ and 

^sO f^a n r^c =: ^o 7^ 0- The overlap domain fto is often referred to as the handshake or 

l^ blending region. In an AtC method, we use the atomistic description on fia and the 

continuum description on Qc- The problem then is how to couple the two different 
descriptions of the material over the overlap region, flo- 
"^ Attempts at this problem thus far can be characterized as either energy-based, 

yr^ where a coupled energy is defined to be minimized; or force-based, where internal 

I and external forces in fij, and Q^ a-re equilibrated. In either case, the resulting AtC 

:• methods often involve some form of blending of the energy and/or the forces over the 

. ^ overlap region. 

j^ The extension of the Arlequin method [2] and quasicontinuum method [13, 17] 

J-j are examples of blended energy AtC methods in which the continuum and atomistic 

energies are combined over fig using a partition of unity. The blended functional is 

then minimized over fl subject to a constraint expressing equality (in a suitable sense) 
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of the atomistic and continuum displacements in flo- 

A standard way to define an AtC method using force blending is to start from the 
variational form of the atomistic and continuum models and blend the corresponding 
weak forms over ilo- We refer to [1] and [14] for investigations of blended force-based 
AtC formulations. A more extensive investigation of these and other AtC methods 
has been presented in [18]. 

In this paper, we formulate and analyze an optimization-based AtC method which 
differs significantly from the existing AtC approaches. The main idea is to cast AtC 
as a constrained optimization problem with virtual Dirichlet controls on the interfaces 
between the atomistic and continuum subdomains. The objective in this optimization 
problem is to minimize a suitable norm of the difference between the atomistic and 
continuum displacement fields over the overlap region Qq, while the atomistic and 
continuum force balance equations on fla and fl^ provide the constraints. Because we 
consider splitting the original problem over 17 into two subproblems over fla and flc, 
we also must impose some form of boundary condition on the two interfaces between 
ria and flc in order to have a well-posed problem. In the context of the optimization 
formulation, these boundary conditions act as Dirichlet controls. However, they are 
virtual, or artificial controls, because the boundaries on which they are imposed are 
artificial rather than actual domain boundaries. 

While the AtC methods in [2,7,13,14,17] also involve a constrained optimization 
formulation, they differ fundamentally from the approach developed in this paper. 
Most notably, the former minimize a blended energy functional subject to constraints 
forcing the equality of the atomistic and continuum displacements over fio. In contrast, 
our approach completely separates the two models and minimizes the discrepancy of 
the atomistic and continuum displacements in flo subject to the two models acting 
independently in fla and fie- The reversal of the roles of the constraints and objectives 
in our approach, relative to blending methods, bears some important theoretical and 
computational advantages. It addresses the problem of blending two physical models 
over a shared spatial region by minimizing instead the mismatch of the deformations 
in the overlap region which is less restrictive on the overall formulation. 

The use of optimization and control ideas for AtC further gives rise to a number 
of attractive theoretical and computational properties. For instance, we are able to 
infer key properties of the method from its atomistic and continuum constituencies, 
as illustrated in the proof of Lemma 4.2. This should be contrasted to the force-based 
quasicontinuum operator which fails to be stable in specific norms even though both 
atomistic and continuum force operators are stable [11]. 

The primary computational advantage of the proposed method is that code to im- 
plement the method can be built upon preexisting code for solving individual atomistic 
and continuum problems. Since the core feature is minimizing the difference between 
solutions of an atomistic model and a continuum model, all that is required is a linking 
program between the two algorithms which carries out the optimization. 

Our AtC work follows a number of previous efforts exploring the use of optimiza- 
tion and control ideas for the design of numerical methods [3-5, 15, 16]. Conceptually, 
our approach is closest to the virtual control techniques for heterogeneous domain 
decomposition developed in [12]. In that setting, both domains are modeled using a 
local, continuum (PDE) model, whereas here we are concerned solely with coupling a 
nonlocal atomistic model with a local continuum model. 

Since the main goal of this paper is to demonstrate the application of optimization 
and control ideas to AtC, we formulate and analyze our method using a linearized 
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Lennard- Jones type atomistic model [8,10,19]. For completeness, we review this 
model in Section 2. We present the new method in Section 3 and analyze its error in 
Section 4. Section 5 summarizes our conclusions. 

2. Preliminaries. This section establishes the notation and defines the model 
atomistic problem. Application of the Cauchy-Born rule [10] to this problem yields 
the continuum formulation. 

2.1. The atomistic model. We consider a chain of A^ + 1 atoms with reference 
(undeformed) positions Xi, i = 0,...,N. The atomic positions in the deformed 
configuration are Xi, i — 0, . . . , N , and Ui = Xi — Xi is the displacement of atom i. We 
assume each atom interacts with its first and second neighbors through a linearized 
Lennard- Jones type potential. Thus, we can effectively think of atoms interacting 
with first and second neighbors via linear springs with spring constants fci > and 
k2, respectively, and equilibrium lengths i and 2i, respectively. For the linearization 
of typical interatomic potentials such as Lennard- Jones, the second neighbor spring 
constant satisfies k2 < [10] and so is not a physical spring, but we will assume 
A:2 < in the following. We also assume nearest neighbor interactions dominate 
second neighbor interactions with the hypothesis 

fci+4fc2>0, (2.1) 

which is also necessary and sufficient for the stability of the atomistic problem [11] 
we consider below. For simplicity, we set the lattice parameter £ = I. 

Under these assumptions, the computational domain is Q := [0, N] D Z, and the 
lattice displacements u — {uq, . . . , un} are elements of the space 

with inner product (•, ■)i2m\ and norm || • H^a/jj) — (•, ■)i2(q)- The left and right "bound- 
aries" of fi are F" = {0, 1} and F+ = {N-1,N}, respectively\ and n° = [2,iV-2]nZ 
is the interior. The size of a domain is ]•], for example, \^\ = (A^-l-l) and \il°\ — N — 3. 
The potential energy of the lattice is the sum of first and second neighbor inter- 
actions 

^"^ k ^"^ k 

i=0 i=l 

We impose homogeneous Dirichlet boundary conditions by fixing the atoms in F~ and 
r+ . The corresponding homogeneous space of admissible displacements is then 

Z^o := {m e W : M = on F- U F+} . 

We assume that a dead load external force, / e Wg, is applied at each atom site 
resulting in a total energy of 

C*H = f"N-(./,"),2(f,). (2.3) 

An equilibrium configuration of the lattice under the dead load is then given by 

m'' = argmin£*°*(u). (2.4) 

u£Uo 



^Thc reason for fixing two boundary atoms is to ensure all unconstrained atoms have a full set 
of neighbors to intcraet with and avoid boundary defects 



The Euler-Lagrange equations 

dui 



0, icn°, (2.5) 



for (2.4) give the force balance constraints at each internal atom. We express these 
constraints using the finite difference operators Ai, A2 :IAq ^^ Uq defined by 

(Aim). = Mj_i - 2Mi + Mi+i, i e fi°, 
(A2m)j = u,j-2 - 2ui + Ui+2, i e ^°- 

From (2.5), the internal force at site i for i € ri° equals {ki/S.iu -\- k2^2u)i. Thus, the 
necessary conditions for the equilibrium of the atomistic system are 

-(fciAiw^ + fcaAaw"), = /», « G 0°, (2.6) 

ui =0, ier-ur+. (2.7) 

The system of linear algebraic equations (2.6)-(2.7) represents the fully atomistic 
problem, which we write compactly as: 

find ii" e Uq such that Am" = / , (2.8) 

where A :— — fciAi — k2^2- 

2.2. The continuum model. To derive the continuum (local) model, we use 
the Cauchy-Born rule ui « l/2('Ui_|_i + Ui_i); see [18], to approximate the second 
neighbor interactions by first neighbor interactions 

iui+i - Ui^iY « 2(ui+i - Ui)'^ + 2{ui - Ui^{f. (2.9) 

Substitution of the Cauchy-Born approximation (2.9) into the atomistic energy (2.2) 
yields the continuum potential energy 



£"= = - ^ fcc(wi+i - Uif - k2(ui - uof - k2{uN - Mw-i)^ , (2.10) 

4=0 

where kc = ki + Ak2 ■ We note that a surface Cauchy-Born correction is not needed 
for (2.10) since we are assuming that Ui = for i G F" U F". 

We now define the total continuum energy under a force / G IAq as 

£*"*(«) -^'(«)-(/,«).2(o)- (2.11) 

An equilibrium configuration of the continuum model minimizes the total energy: 

u'= = argmin£*°*(u). (2.12) 

The Euler-Lagrange equations for (2.12) are given by 

-(fc.Aiu^), = /„ zGl}°, (2.13) 

u1= 0, iGF"UF+. (2.14) 

The system of linear algebraic equations (2.13)-(2.14) represents the continuum prob- 
lem. Setting C = — fc^Ai, this system assumes the form: 

find u" G U(> such that Cu" = f . (2.15) 
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2.3. The continuum modeling error. The error in the approximation of the 
atomistic solution by the continuum solution is given by the following proposition. 

Proposition 2.1. There exists a fixed constant co, independent of N and u°', 
such that 

Wu"" - u^Wpi^) < coN^WAfu'^Wp^n). (2.16) 

Proof. To estimate the continuum modeling error, we observe that 

A-C = - (fciAi + fcaAa) + (fci + 4fc2)Ai = -fcalAa - 4Ai) = -fcsA?. 

Now C is just the ID discrete Laplacian on Q with homogeneous Dirichlet boundary 
conditions at atom sites 1 and iV — 1. So, the minimum eigenvalue for C is Ai = 
4fcc sin I „, '^ ., ) where n = iV — 3 is the number of unconstrained atoms. Using that 
u°- —u'^ = at atoms 1 and A^ — 1 implies Cu^ = / = Au"-, which yields the following 
bound: 



ulmn)< \\C-^Pin)-\\Ciu''-n')\\e^n) 






Akr sin^ 



2(n+l) 



D 



3. Optimization-based AtC formulation. As with any AtC formulation, we 
begin by splitting 17 into atomistic, continuum, and overlap regions 

i7a = [o,i]nz, rtc^[K,N]r\Z, no^^ar\nc = [K,L]nz, 

where < K < L < N. The strict interiors of these domains are 

17° = [2,L-2]nz, n° = [K + 2,N-2]nz, ni = ninn° = [K + 2,L-2]nz, 

and their boundaries are 

r- = {0,l} and r+ = {i-l,L}, 
r- ={K,K + 1} and Tf = {N -1,N}, 
T- = {K,K + 1} and r+ = {L-l,L}. 

Here and in the remainder of the paper, we find it convenient to use modified Vino- 
gradov notation where the implied constant is independent of the parameters K, L, 
and N. Thus, X > Y means there is a positive constant c such that X > cY with c 
independent oi K,L, and N. 

Recall that the main objective of an AtC method is a stable, accurate, and ef- 
ficient approximation of the lattice displacements by using the atomistic model on 
r^a, employing the continuum approximation on Clc, and accurately merging them to- 
gether on flo- It follows that the efficiency of AtC methods hinges on the assumption 
that the atomistic region is small compared to the continuum region. On the other 
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hand, it is intuitively clear that a stable and accurate AtC method requires some 
conditions on the size of JIq. These assumptions are pivotal to our analysis and we 
formalize them below. 

Assumption A. There exists a real number p > 1 such that 

L < N^/P. 



Assumption B. There exists a real number 7, £ < 7 < 1, such that 

L-K 

= 7- 



L 

In other words, we assume that \^a\ ^ |0|^/p and \Vlo\ — 7|f^a| is such that |51o| > 3, 
i.e., the overlap region's size is at least twice the size of the interaction range; see 
Fig. 3.1. The assumption that 7 is constant means the ratio of the overlap width to 
the size of the atomistic region is constant, or cquivalently, that the ratio oi K to L 
is constant. 



^n 



n„ 



fi. 



K 



L 



N 



Fig. 3.1. Decomposition of Q. Squares are in the atomistic region Qa, circles are in the 
continuum region Qc, and crosses are in the overlap region ilo- 

We formulate the optimization-based AtC method in two steps. The first step 
defines independent atomistic and continuum subproblems on ila and Qc, respectively, 
whereas the second step merges these problems by minimizing the mismatch of their 
solutions on flo- To describe the first step, we introduce the spaces 



Ua 


= {u 


Oa- 


-> R 


u- 


= Oonr-}, 


l^afi 


= {u 


^a- 


-> M 


u- 


= Oonr-ur+}, 


We 


= {u 


^c- 


->M 


u = 


= on r+} , 


Uc.^ 


= {u 


^c- 


->M 


u = 


= on r- u r+} , 



for the subdomain displacements and the "trace" spaces 



Aa = {. 



i^ and Ac ~ {f 







(3.1) 



for the displacement values on the artificial domain boundaries r+ and F^ . We denote 
the standard l"^ inner product and norm on these spaces by (•, ■)t,2i„\ and || • ||^2(ct), 
where cr stands for the appropriate domain under consideration. The trace spaces 
provide the boundary conditions on F+ and F^ necessary to formulate well-posed 
atomistic and continuum problems on $!„ and l^c- 

Let Aa and /° be the restrictions of A and / to the interior il°. Likewise, let Cc 
and f^ denote the restrictions of C and / to the interior f7° . The local nature of the 
continuum subdomain operator Cc necessitates the need for only a single boundary 
constraint on T~ whereas the atomistic subdomain operator Aa requires two con- 
straints on F+. For this reason, when referring to the continuum model, we adjust 
the definitions of f^c, ^°i T^, and F+ to be 



VLr 



[x,Af-i]nz, f7° = [x + i,Ar-2]nz, 



{K}, Tt = {N-l}, 



with analogous changes made to overlap boundaries and interiors and the displacement 
and trace spaces. Thus, the left, artificial boundary and the right, true boundary of 
the continuum region are single atoms. The continuum displacement at atom N is 
zero since this is a true boundary condition of the original problem. To simplify 
notation in upcoming computations, we set N :— N ~ 1; see (3.2). 



rr rj r+ r 

N N 



V • • • • • X O O ox X ••••••XX 

1 , K , L — 1 L 



n° 



Fig. 3.2. Trace spaces and interiors oJQ,a, ^c- The interior of Qo is depicted with open circles. 

We define the atomistic subproblem as a restriction of (2.8) to fia with inhomoge- 
neous boundary conditions at the artificial atomistic boundary r+, i.e., given 6°- g A° 
we seek u°- E Ua such that 

f Aau" ^ f on n° 

1 ! • (3.2) 

y u" = e" onT+ 

Similarly, the continuum subproblem is a restriction of (2.15) to f^c with inhomoge- 
neous boundary condition at the artificial continuum boundary F^ : given 0^ G Ac we 
seek M^ e Uc such that 

r Ccu" = f on ni 

1 • (3.3) 

[ m'^ == 6*'= on r- 

Thanks to the boundary conditions prescribed on the artificial boundaries, the sub- 
domain problems (3.2)-(3.3) are well-posed and can be solved for any given 0" and 
9'^. However, because 6*" and 9'^ are unknown, the solutions to (3.2) and (3.3) cannot 
yet be determined. 

The second step in the formulation of our AtC method is the merging of (3.2) 
and (3.3) into a single well-posed problem for the unknown states w" and u'^, and the 
unknown boundary conditions 0" and 9'^. Intuitively, we desire that 

r^u'^onr-, 0"wu'=onr+, and u" « u'^ in f7°. (3.4) 

In many hybrid AtC methods these, or similar conditions, are used to constrain the 
hybrid force balance equations or the minimization of a hybrid energy functional; see 
e.g., [2,7, 14]. However, there is no canonical way of enforcing strong or weak type 
equality of fundamentally different atomistic and continuum solution states. 

The cornerstone of our optimization-based AtC approach is to view (3.4) as the 
optimization objective rather than as the constraint. Specifically, in the context of 
our model problem, the quantity 

II"" - ^%Hn.) - W^"- m^A.) + II"" - ^%Hni) + U - ^1I'^(aj (3.5) 

provides a notion of an artificial "mismatch" energy between the solutions of (3.2) and 
(3.3) in the overlap region. Instead of forcing this energy to be exactly zero, which does 

7 



not yield a problem with a solution, we seek to minimize it subject to the atomistic 
and continuum force balance equations (3.2) and (3.3) holding independently in fla 
and flc- Succinctly, our new AtC formulation is the following constrained optimization 
problem: 

1„ „ .,,2 (AaU^^f on 17° C^u'^^f on 17° 



min -\\u°- - u'^WJ ,„ ^ s.t. < 
.,„se»,e=}2" "^^(^^°) I u« = 0« onr+ 



{«»,«-,«» ,e<=} 2" ""^^"°^ ■ ' 1^ u«==6'" onr+ ' 71^= == ^ on T" 

In the language of constrained optimization, the functions u" G Via and u'^ G Z^c are 
the states, and the artificial boundary conditions 0" G Ka and 0^ G Ac are the controls. 
The purpose of the controls is to allow the states to adjust so as to provide the smallest 
possible value of the objective while still satisfying the constraints. In the context of 
(3.6), 0" and 6'^ are virtual boundary controls, as the boundaries r+ and F^ are an 
artifact of the domain decomposition into atomistic and continuum parts. 

We shall show below that the optimization problem (3.6) is well-posed. But before 
investigating this, we show the optimization-based AtC formulation (3.6) satisfies a 
patch test criterion. 

3.1. Patch Test Consistency. The bane of all atomistic-to-continuum coupling 
mechanisms is the existence of nonphysical ghost forces arising on the interface of the 
contimmm and atomistic regions [9,20]. The patch test is a well-known test for 
determining the existence of ghost forces by checking whether a uniform strain is an 
equilibrium solution to the proposed method on a perfect lattice under zero external 
forces [18,20]. As with force-based methods, the optimization formulation (3.6) is 
patch test consistent by design. Indeed, if we replace the homogeneous Dirichlet 
boundary conditions by the inhomogeneous boundary conditions 

ug = 0, ul^F and u%_^ ^ {N - l)F, u% = NF, 

where -F > defines a macroscopic displacement gradient, and if we take take /" = 
0, f = 0, then it is straightforward to verify that (3.6) has a minimum of achieved 
when u^ = iF and u1 — iF . This is due to the fact that both atomistic and continuum 
operators are patch test consistent individually. 

3.2. Well-Posedness. To establish that the optimization-based AtC formula- 
tion (3.6) is well-posed, we switch to the reduced space form of the optimization prob- 
lem, which requires the elimination of the states from (3.6). In our case, this task is 
trivial because for any pair of virtual controls {6°',0'^} G A^ x A^ the constraints 

: f on 17°, r Ccu'^ ^ r on f7°, 

and J " •' -' (3.7) 

:r onr+, [ u'^^e'' onT^, 

have unique solutions u"" = m" [O"") G Ua and u^ — u^ {6^) G Uc- Using these solutions 
in (3.6) transforms the latter into an equivalent unconstrained minimization problem 

min \\\u''{e'')-u-{e-)\\%(^y (3.8) 

This problem, in terms of the virtual controls only, is the reduced space form of (3.6). 
We analyze (3.8) following the strategy in Gervasio et al. [12]. Specifically, for any 
given {0°- , 6'^} G A^ x A^ wc split the solutions of the constraint equations (3.7) as 

m"(6'") = v^iO") + m"'° and u"(0") = v^iO") + u"'° (3.9) 




where the homogeneous components u"'" e Ua.o and u'^'" e Uc,o solve 

Aau'''° = /" and Ccu"^" = /', (3.10) 

respectively, whereas v°'{6"-) E Ua and v'^{6'^) e Wc solve 

f Aa-y" = on r2° f Ccw'^ = on f7° 

<^ " and <^ " , (3.11) 

[ v" = 6" on r+ [ i;'^ = 9" on T^ 

respectively. We will prove the following stability result for (3.11) in Appendix A. 

Lemma 3.1. For any {0°-, 9"} G A^ x A^ the solutions v°-{0°-) and v''{9'') to 
(3.11) satisfy the bounds 

h%n\\%in.) < L\r\\%(r+y 

(3 12) 

Remark 3.1. The bounds (3.12) continue to hold when we take homogeneous 
Dirichlet boundary conditions on F^ and T~ in (3.11) and inhomogeneous boundary 
conditions on F^ and F^t by replacing F J with T~ and F~ with F^t . 

Using the decomposition (3.9), the reduced space problem (3.8) assumes the form 

,, min J||«"(e")-i'^(mil^(o„) + («"r)-«'r),""'"-"''°),.(a) 

{e»,e<^}eA„xA, 2 ^ °> /« is«oj /„ -|o\ 

The following result is key to proving that the reduced space problem (3.8), respec- 
tively (3.13), has a unique minimizer. 
Theorem 3.2. The form 

{{e\e-},{^i\^,^}) :=. (^;'^(r)-^;^(n,«'^(A*'^)-^'^(M^)),.(n„) (3.14) 

defines an inner product on A^ x Ac. 

Proof. The proof follows from the inequality (4.15) in Lemma 4.2. 

D 

Theorem 3.2 allows us to recast the reduced space problem (3.13) as 
1, 



1, 



(3.15) 



\ afi „c,0||2 



where || • |J£*(a„xAc) i^ the norm induced by (3.14). The necessary optimality condi- 
tion (Euler-Lagrange equation) for (3.15) is the following variational equation: find 
{6*", 6"} e Aa X Ac such that 

{{e\e^},{^i\^,^}) = - {u--" -u^'^,v^{p'^)-v%^i'^))^,^^^^^ (3.16) 

for all {ij,"-,^'^} e (Ad X Ac) . Theorem 3.2 and the Riesz representation theorem im- 
ply that (3.16) has a unique solution, thereby establishing the well-posedness of the 
reduced space problem (3.13). 
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Remark 3.3. There are two principal pathways for the analysis and the numerical 
solution of the constrained optimization problem (3.6). The first one, which we adopt 
in this paper, relies on the strictly convex reduced space problem (3.8) and its equivalent 
forms (3.13) and (3.15). In this case the practical implementation of the AtC method 
involves the solution of the strongly coercive Euler- Lagrange equation (3.16), followed 
by the recovery of the atomistic and continuum .states from the virtual controls. 

The second pathway relies on Lagrange multipliers to enforce the constraints in 
(3. 6) and yields a saddle-point optimization problem. The latter can be analyzed using 
the Brezzi's theory [6], while implementation of the AtC method then requires the 
solution of a weakly -coercive, mixed-type optimality system. We note that showing the 
conditions of the Brezzi theory is essentially equivalent to showing the well-posedness 
of the reduced space problem. 

Remark 3.4. Let A^ be the restriction of A to the continuum subdomain 51^. 
Substituting Ac for Cc in (3.6) yields a constrained optimization formulation that 
is equivalent to the global atomistic problem (2.8). The corresponding reduced space 
problem and its Euler-Lagrange equation differ from (3.15) and (3.16) only by the use 
of the atomistic operator Ac instead ofC'c. As a result, the variational problem (3.16) 
can be thought of as resulting from a modification of the bilinear form associated to 
the original problem very similar to a nonconforming finite element method. 

4. Consistency and Error Analysis. This section analyzes the error be- 
tween the true solution u" to the atomistic problem (2.6) and the solution of the 
optimization-based AtC formulation (3.6). 

To proceed with our analysis, let {0°, ^^p} e Aq x Ac denote the optimal solution 
of the reduced space problem (3.15) or, what is the same — the solution of the Euler- 
Lagrange equation (3.16). The optimal solution of the full problem (3.6) is then given 
by {Kp^ <p: ^op> Qlp}^ where 

<p = ^'"(^op)+""'° and <p = «^(0^p) + u^^o. 

These are the optimal states. Using these states we define the AtC approximation to 
w" as 

yulp mVLc\^o- 
For the error analysis, it is convenient to express the approximate AtC solution as 

where the affine operator, P : Aq x Ac ^ Wq, is defined by 

P({/x«,A^=}):-| ""'° + «"(M")in^a, V^, A^^} G A, x Ac . (4.2) 

Thus, the error of the AtC approximation (4.1) is 

\\u^ - «"*1l£^(o) ^\\u--p {{e:^, ei^}) |i,.(j,) . (4.3) 

To analyze (4.3) it is advantageous to split P into a linear part, Q, and constant 
term, U^ , dependent only on the homogeneous data, i.e., P = Q + U^ where 



v'^(m") in r!„, 

w=(/i^) in r2c\^o, ~ " ' I «'=■" in r2c\^o . 



Q(K,M^}):= Tc) c. : '. and [/' .- , ^„ 
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We also introduce the trace operator r : Wo i-^ Aa x Ac such that 

rH:={("^^'),M}:-KH,r^(")} Vn G Wo . (4.4) 

Because r"(u") contains the exact values of the atomistic solution, it follows that 

«1n„-«"(r"(""))+w"'°- (4.5) 

We define the continuum lifting of the exact atomistic trace on Aj, as 

u'=:=i;<=(r'=(u")) + M^'0. (4.6) 

It is a matter of unraveling these definitions to see that 

{v'^ in O 
c • nln ■ (4.7) 

u m sic\"o 

To estimate the AtC approximation error we split (4.3) into two parts: 

= \\u^-P {r{u^)) + P {r{u^)) P {{ei^, ei^]) ||,.(o) 

< 11^" - P(r(^"))||,.(o) + \\Qr{u'^) + U'-Q {e^,p,e^,^} - U"\U.^n) (4.8) 
= 11^" - P(r(^'^))||,.(o) + \\Q {r{u'^) - {O:p,e'op}) ll^^(n) 

< 11^" - P(r(^"))||,.(o) + IIQII • ||r(M«) - {0«p,0^p} ||,*(A„xA.), 

where || • ||^*(a„xAc) is the norm induced by (3.14) and 

IIQII = sup '-^. (4.9) 

{A<».M<^}eA„xAe II iM :M IIU*(AaXA,) 

The first term in (4.8) is the consistency error of the operator P. Using (4.7) 

\\u^ - P (r(u")) ||,2(f,) - ||«'' - u^\U2(n^\n^) , (4.10) 

i.e., the consistency error is confined to the purely continuum region. The second term 
is proportional, up to a factor of ||Q||, to the approximation error^ in the solution of 
the reduced space problem (3.13). We proceed with an estimate of the approximation 
error, followed by a bound on the operator norm ||QII- 

Lemma 4.1. Let u'^ solve (2.6) and {^op,^™} be the minimizer of (3.13). Then 

Hun - {{O^op.O^^p}) ||,*(A„xA.) < ll^t" - u^Wmn^) ■ (4.11) 

Proof. We bound the approximation error directly by noting jS" , ^opj solves the 
Euler-Lagrange equation (3.16) of the reduced space problem. As a result, 

lk(^")-({^op,^op})ll£*(A.xA.) 

{r{u-), {^l-, ^i-}) + («-." - u-.o, v-i^I-) - «^(m^)),.(j,^) (4.12) 



sup 



-T 1 1 r Til 

{m^.m-'I^O II iM''>M''l ll£*(AaXA,) 



^This error measures the difference between traces of the true atomistic solution ii" and the 
approximate AtC solution (4.1). 
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Using Definition (3.14), (4.5) and (4.6), wc can obtain 

= {v^ir'^iu'^)) - v^{r^{un),v^{fin - «^(m^)),.(,j„) 

Using this identity in (4.12) completes the proof. D 

The following lemma estimates the norm of Q. 

Lemma 4.2. Under Assumption A and Assumption B, the norm of Q is bounded 
by 



IIQII<7-Yr^' (4.13) 

where the implied constant is allowed to depend on p. 

Proof. Definition (4.9) implies that (4.13) will follow if we can show that 

||g({Ai",M^})||,.(o)< 7"'y^^llK,A^^}||£MA.xA.) V{M",M^}eA„xA,. 

(4.14) 
On the other hand, the definition of Q and (3.14) imply that (4.14) is equivalent to 

Wv^'il^nWlin.) + lk^(M^)lll^(o./o„) < 1-' {y^) II^"(^'") - «'(A^')ll'^(f^.) (4-15) 

for all {/x",/z^} e Aq X Ac. To prove (4.15) we use the structure oiv"-{fi°-) and ti^(/i^). 
Recall that v'^{ii'^) E Uc solves the continuum submodel 

Cv'^^Q in 17°, v'k^ii'k, v% = Q. 

It is straightforward to verify that ti'^(/i^) is a linear function, i.e., 

N — i 
w-=ac^ , (4.16) 

where ac = nl^. 

On the other hand, u°(/^") G Ua solves the atomistic submodel 

We decompose this field as v°-{ii°') — v"^ + v^ + v^ + v^ where 



f,i=ai- and vj ^ aaA-^^Vi - K (4.17) 

Ij 

are the dominant linear and exponential modes with ai and a2 determined by the 
boundary conditions below. Meanwhile, v^ and v^ are corrections to ensure that 
v'^ifi") = onT", i.e., 

in f}° f Aav"^ =0 in ^° 

(4.18) 
[ w'^ == on r+ 
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inf}° 




r AaV^ = 





in ni 


v^ = 


-v^ 


onr- 


and 


V^^ 


-1.1 


on r, 


v^ = 





onr+ 




V^^ 





on r, 



To obtain v^ we need the roots of the characteristic polynomial of Aa 

p{(t) = -^20-'' - kicr^ + (2fci + 2k2)a^ - ha - fcz, 
which are given by (see [9]) 

fci + 2fc2 ± y^kj +Akik2 



Ai — A2 — 1; A3_4 — 



-2h 



1 



We define w^ by setting A == A4 = fci+2fc2 ^kj+j k^k^^ .^ ^^ _^^^^ ^^^^ ^^^^^ < A < 1, 
as seen from the assumptions fci > 0, fe2 < 0, and fci + 4fc2 > 0. 

The coefficients ai and a2 are uniquely determined from the boundary condition 
ti"(/i°) = /i" on r+, which yields the following 2x2 system: 

1 - 7 ) "1 + AVi - Ka2 - A*2-i, ,,,„, 

■t-/ (4.19) 

q;i + \/L — Ka2 = /i/,. 
Recall that ac — f^^- I^ ^he following, we define 



According to Remark 3.1, the result of Lemma 3.1 applies to v'^ and u^, and so. 

Since w^L- = (0, ai/i) and f^L- = oi2\/L — K (A^, A^^^) , we have the bounds 

\\vXHn.)<>^^^''-'''L{L-K)al and ||«^||,'.(,,„) < i^. (4.20) 

Using (4.20) yields the following lower bound for the right hand side in (4.15): 

II«"(m°) - v^{p")\\t^i-io) >\\v^+v'^ - v^Wpf^no) - \\A\p(i\>) - \\A\p(i-io) 

/ 1 \ (4.21) 

We proceed with estimating 

\\v' +v'- «1l|.(Oo) -\\v'- «1ll(no) + 2 («^ - «^ «'),.(o„) + ll«'lll^(no)- (4-22) 

The term \\v^ — v'^\\'J2(q \ is similar to the term in (3.14) defining the trace norm 
II ' ll£*(AaxAc)i but it is simpler in that both v^ and v"^ solve continuum problems. We 
will prove in Appendix B that 

\\v' ^^111.(00) >iL- K)j' {al + aj) (4.23) 

for large N. Furthermore, summing a finite geometric series shows 

ll« \\p{no) = Y^^ [ ) 2- (4.24) 
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Intuitively, one should suspect the cross term, {v^ — v^,v'^)i2(^fi^), in (4.22) to be 
estimable for large overlap widths since the exponential term u^ is not well approxi- 
mated by any linear function. We now calculate via explicit summation that 



{v'~v-,v^) 



i^iilo) 



y^ ( ai-^ - ac^^ ^ ) a2X^ WL- K 



=K 



L 



N -K 

L 



N~i , 
N -K' 



. L—i 



= a2\/L - K ai ^ -A^ ' - etc Y^ 

V i=K i=K 

L 

> - \a1a2\VL -kY^ x^-' - \aca2\VL -kJ2^ 



(4.25) 



L-i 



i=K 



i=K 



> -Vl-k 



1 - A^-^+i 
1-A 



Using the inequalities (4.23), (4.24), and (4.25) in (4.22) produces 
\\v'+v'-vlUno)^{L-Kh'{al+al) 



+ {L-K) 
(L-K) 



I _ y2{L-K+l) 



Ot — Vi — K 



1-A2 ,'"2 

1 1 - A^-^+1 



1 - A^-^+i 



VL-K 1-A 

1-A2 



1-A 

("c +ai) 
1 1 - A^-^+i 



(4.26) 



^/L-K 1-A 

For a sufficiently large overlap region, there holds 

r- ri-A^-^+1 1 + A 1 

which guarantees the positivity of the terms multiplying a^ + a^ and a^ above. Then 
since 7^ < 1, we obtain from (4.26) that 

\\v'+v'-v%.^,,^^>{L-K)j^a\ (4.27) 

Similarly, using (4.27) in (4.21) yields 

\\v''-v'\\p^n„)> VL^^ja. (4.28) 

To complete the proof, we use the above results to estimate the left-hand side in 
(4.15): 

I _ y2(l + L-K) 



<A[\\v 



<A{La{ + {L-K)- 

+ {N-L)al 



1-A2 



+ L{L-K)X'''^-^al + ^ 
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Estimates of the norms of v^, v"^, v^, v , and v'^ follow from (3.12) in Lemma 3.1, 
(4.24) and (4.20), respectively. The final inequality, which establishes the assertion of 
the lemma, is a consequence of (4.28). D 

All results necessary for the completion of the AtC approximation error bound in 
(4.8) are now in place. 

Proposition 4.1. Let u"- solve (2.6) and {Ogp.d^p} be the minimizer of (3.13). 
The AtC solution satisfies the error bound 



ll^"-«''*1l£^(o)< fl + 7-y^^j Wu^-W^Wmn.y (4.29) 

Proof. Recall the split of the AtC solution error into a consistency error due to 
P and the approximation error in the reduced space problem (3.13): 

Wu^'-u'^'lpm - 11^" - p {{e^p, 0Sp}) Wi^n) 

< \\u^- P{r{u^))\\,.^n) + IIQII • Mun - {O^op.O^p} II£*(a„xA.), 

Using (4.10) for the consistency error, (4.11) for the approximation error, and (4.13) 
for operator norm yields the result of the proposition. D 

Proposition 4.1 reveals that the accuracy of the AtC approximation is determined 
by two independent factors. Replacing the atomistic model with a continuum model 
on flc introduces the continuum modeling error ||u" — u'^\\p(^Q^-), which is independent 
of the choice of the coupling mechanism. An inherent assumption in atomistic-to- 
continuum coupling is that the continuum model closely approximates the atomistic 
model in the continuum region, which can be expected when there are no defects in 
the continuum region [18]. Thus, we expect \\u°' — u'^|!£2(f2e) ^o be small so long as this 
assumption holds. On the other hand, the coupling mechanism via the optimization 
framework introduces the prefactor 




which depends on the size of the overlap region. As can be expected, the AtC error 
is inversely proportional to the size of flo- 

We can precisely estimate the modeling error by applying the estimate (2.16) 
to the domain flc- The operator C is now the ID discrete Laplacian on flc with 
homogeneous Dirichlet boundary conditions at K and iV — 1. Thus, the minimum 
eigenvalue for C is Ai = Akc sin^ I 2(n+i) ) where now n^N — K — 2 is the dimension 
of C. We then have since u° — u'^ — at atoms K and iV — 1 that 

ll«"-"1l£^(o.)< {N ~ Kf\\Alu-\\,2^^^y (4.30) 

The estimate (4.30) confirms that the modeling error is small whenever u'^ is smooth 
over the continuum region in the sense that ||A^u"||^2(f2^) is small. 

By using the modeling error bound (4.30) in (4.29), we obtain the following the- 
orem for the AtC error estimate. 

Theorem 4.2. Under the conditions of Assumptions A and B, let u" solve (2.6), 
and let l^™,^™} ^^ ^^^ minimizer of (3.13). Then 



u'''^hHn)< (l + 7-S/;^j(^-^)'l|A?Sl|,.(o.), (4.31) 
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where u"'^'^ — P{{(^opt(^op\) ^^ ^^^ ^^C! solution and P is defined by (4.2). 

Note that the dependence of the AtC error on the size of the overlap domain 
is unchanged, i.e., as the overlap width is increased, the error decreases. In the 
present situation, we did not coarse-grain the continuum region, so the only increase 
in complexity comes from increasing the size of the atomistic and continuum regions. 

Thermodynamic limit. By letting A^ -^ oo, the problem above is an example of 
a thermodynamic limit. A further estimate would typically be obtained by assuming 
the fully atomistic solution decays sufficiently rapidly as A — )■ oo. Sec [18] for an 
analysis in this setting for quasicontinuum methods. 

We may conversely introduce an interatomic spacing parameter e with e dependent 
norm 



Setting e — N^^ maintains 51 == [0, 1] and scaling the lattice by j h^ ej scales 

ll^iw||i^2(fj^) H^ e''||Aiu||^2(j2^)^,. 
The estimate in (4.31) under this scaling is 

3 

\\u''-u'''lpm,e< ^^L=||A?^'^||,.(^^),„ (4.32) 

which is the scaling limit as e ^ 0. See [11] for the derivation of a similar estimate 
in the case of the force based quasicontinuum operator in which e is maintained as a 
parameter throughout. Recalling Assumptions A and B, if the overlap region f2o has 
width \flo\ := [L — K)e^/P in the scaling limit, then we obtain the bound 

- + — 
||7i"-z.«*=||,.(f,),,< "l^W^lu-Wp^n.y,. (4.33) 

|S2o|2 

Hence, we may achieve any power of e in the interval (|, 2j for p > 1. 

5. Conclusion. This paper formulates and analyzes a new, optimization-based 
strategy for atomistic-to-continuum coupling. Specifically, we pose the problem of 
coupling a non-local, atomistic description of a material with a local, continuous 
description as a constrained optimization problem. The objective is to minimize the 
£^ difference between the continuum and atomistic displacement fields over an overlap 
region, subject to constraints expressing the atomistic and continuum force balances 
in the respective subregions. The traces of the atomistic and continuum solution 
components on the boundary of the overlap region act as virtual boundary controls. 
Thus, our approach can be viewed as an extension of the heterogeneous decomposition 
method [12] to the AtC context. 
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Appendix A. Stability of atomistic and continuum problems. In this 
appendix, we prove the result stated in Lemma 3.1: 



,,a /'^a\ II 2 <^ T \\ 



\\v^{e^)\\ 



90112 
Ilf2(r+)' 

3C||2 



(A.l) 



The second bound is a direct consequence of the maximum principle for the continuum 
operator C — — fccAi. Recalling that v'^{9'^) is zero on r+ and equal to O^^ on T~ = 
{-fC}, we have 



N-l 



N-l 



\v'mf 



e.^(na) 



E «)^ ^ E 



'K) 



{N - K){e\ 



K) 



(A.2) 



=K 



^K 



To prove the first bound in (A.l), we note that the atomistic solution, u°(6'"), 
may be written as 



^'"(^")n = A ^ + /32^^ + /53A" + /34A^-" 



(A.3) 



where < A < 1 was defined in Lemma 4.2 and the coefhccnts arc determined via the 
boundary conditions v°-{9"-) = on F^ and v°-{9°-) — 9°- on r+. Specifically, 



Tl 



//3A / 

/32 



i_ 

L 

V 1 



1 

L-1 

L 
J_ 

L 




1 



^^^ 



A A^-i 
A^-i A 

A-^ 1 y 



(h\ ( \ 



/32 
/33 

V/^4y 



L-l 

v^2 y 



(A.4) 



where 



/O 1 1 0\ 

1 A 

1 A 
yi 1/ 



=: T as L — > 00. 



(A.5) 



For any (5 > 0, we can therefore choose L such that ||T^ "'^H < \\T '^W +5, and hence 

{Pi + /3| + /3| + /3|) < (r-^ll + 6f ({9l_,f + (0£)^) . (A.6) 

Using successive Cauchy inequalities and explicit summation of finite geometric series 
yields 

ii«''r)ii£^(oj 



< 4 PiL + P'^L + /33^ 



1-A 



L+l 



1-A 



-fit 



A-i-A^ 
A-i-1 
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for large enough L by (A. 6). 

Appendix B. Estimate of \\v^ — u^||^2(f2j,). Finally, we establish the estimate 
(4.23) under the conditions of Assumptions A and B. Recall that v'^ and v^ are defined 
in (4.f6) and (4.17), respectively, and so, 

^ / N — i z \^ 

^'-^''ll'^(o„) = I] (ac^^3-^-ai-j ^Aal~2Ca,ai + Bal (B.l) 

i—K ^ ^ 

where the coefficients of the quadratic form in (B.l) are given by 

respectively. Summing the finite series for each coefficient yields A — (3 ■ A, B ~ 13 ■ B, 
and C = P ■ C where the common factor is/3= (1 + L — K) and 

6iV2 + 2^2 + 2K^ - 6KN - 6LN + 2KL + L-K 

2L^ + 2K^ + 2KL + L-K 
B = 



C = 



6L2 
2^2 + 2K^ + 2KL - 3KN - 3LN + L-K 



6L{K - N) 

respectively. Using that K — {1 — j)L allows us to further write the coefficients as 

N^ + L^l-j+ 17^) + i7£ - LN{2 - 7) 
(iV-(l-7)i)2 

LO^T+M+h iv(i-h)-h-^(i-7 + h^) 

L ' 7V-(1-7)L 

Assumption A implies that limjv_5.co L/N = 0, and so 

A^Aoo^l, B^Boo = l-7+57^ and C ^ Coo ^ I - h ■ 

Let < Ai < A2 be the eigenvalues of the quadratic form AooCXc ~ 2CooQ!cai + BaoOc{. 
Using the expressions for the determinant and the trace of the quadratic form, A1A2 = 
AooBoo ~ C*^ and Ai + A2 = A^o + Boo, we can estimate 

A1A2 ^ A1A2 ^ kL__ > 1^2 

A2 - A1 + A2 1 + (1-7+172) - 24' 

This completes the proof. 



18 



REFERENCES 

[1] S. Badia, M. Parks, P. Bochcv, M. Gunzburgcr, and R. Lehoucq. On atomistic-to-continuum 

coupling by blending. Multiscale Modeling & Simulation, 7(1):381— 406, 2008. 
[2] P. Baunian, H. Ben Dhia, N. Elkhodja, J. Oden, and S. Prudhomnie. On the application of 
the Arlequin method to the coupling of particle and continuum models. Computational 
Mechanics, 42:511-530, 2008. 10.1007/s00466-008-0291-l. 
[3] P. Bochev and D. Ridzal. Additive operator decomposition and optimization— based reconnec- 
tion with applications. In I. Lirkov, S. Margenov, and J. Wasniewski, editors. Proceedings 
of LSSC 2009, volume 5910 of Springer Lecture Notes in Computer Science, 2009. 
[4] P. Bochcv and D. Ridzal. An optimization-based approach for the design of PDE solution 

algorithms. SIAM Journal on Numerical Analysis, 47(5):3938-3955, 2009. 
[5] P. Bochev, D. Ridzal, and D. Young. Optimization-based modeling with applications to trans- 
port. Part 1. Abstract formulation. In I. Lirkov, S. Margenov, and J. Wasniewski, editors, 
Proceedings of LSSC 2011, Springer Lecture Notes in Computer Science, Submitted 2011. 
[6] F. Brczzi. On the existence, uniqueness and approximation of saddle-point problems arising 

from lagrangian multipliers. RAIRO Anal. Numer, 8(2):129-151, 1974. 
[7] L. Chamoin, S. Prudhomme, H. Ben Dhia, and J. Oden. Ghost forces and spurious effects in 
atomic-to-continuum coupling methods by the Arlequin approach. International Journal 
for Numerical Methods in Engineering, 83(8-9):1081-1113, 2010. 
[8] W. Curtin and R. Miller. Atomistic/continuum coupling in computational materials science. 

Modelling Simul. Mater. Sci. Eng., 11:R33-R68, 2003. 
[9] M. Dobson and M. Luskin. An analysis of the effect of ghost force oscillation on quasicontinuum 
error. Mathematical Modelling and Numerical Analysis, 43:591-604, 2009. 

[10] M. Dobson, M. Luskin, and C. Ortner. Accuracy of quasicontinuum approximations near 
instabilities. Journal of the Mechanics and Physics of Solids, 58:1741-1757, 2010. 

[11] M. Dobson, M. Luskin, and C. Ortner. Stability, instability, and error of the force-based 
quasicontinuum approximation. Archive for Rational Mechanics and Analysis, 197(1):179- 
202, 2010. 

[12] P. Gervasio, J.L. Lions, and A. Quarteroni. Heterogeneous coupling by virtual control methods. 
Numerische Mathematik, 90(2):241-264, 2001. 

[13] B. Van Koten and M. Luskin. Analysis of energy-based blended quasi-continuum approxima- 
tions. SIAM J. Numer. Anal, 49(5):2182-2209, 2011. 

[14] X. Li, M. Luskin, and C. Ortner. Positive-definiteness of the blended force-based quasicontin- 
uum method. SIAM J. Multiscale Modeling & Simulation, 10, 2012. arXiv:1112.2528vl. 

[15] J.L. Lions. Virtual and effective control for distributed systems and decomposition of every- 
thing. Journal d'Analyse Mathematique, 80:257-297, 2000. 10.1007/BF02791538. 

[16] J.L. Lions and O. Pironneau. Virtual control, replicas and decomposition of operators. C. R. 
Acad. Sci. Paris, 330(l):47-54, 2000. 

[17] M. Luskin, C. Ortner, and B. Van Koten. Formulation and optimization of the energy-based 
blended quasicontinuum method. Computer Methods in Applied Mechanics and Engineer- 
ing, 253:160-168, 2013. arXiv: 1112.2377. 

[18] Mitchell Luskin and Christoph Ortner. Atomistic-to-continuum coupling. Acta Numerica, 
22:397-508, 2013. 

[19] M. Parks, P. Bochev, and R. Lehoucq. Connecting atomistic-to-continuum coupling and domain 
decomposition. SIAM J. Multiscale Model. Simul, 7(l):362-380, 2008. 

[20] V. B. Shenoy, R. Miller, E. B. Tadmor, D. Rodney, R. Phillips, and M. Ortiz. An adaptive 
finite element approach to atomic-scale mechanics-the quasicontinuum method. J. Mech. 
Phys. Solids, 47(3):611-642, 1999. 



19 



